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Abstract 

We present a scheme for the calculation of linear optical properties by the all-electron full- 
potential linearized augmented planewave (LAPW) method. A summary of the theoretical back- 
ground for the derivation of the dielectric tensor within the random-phase approximation is pro- 
vided together with symmetry considerations and the relation between the optical constants. The 
momentum matrix elements are evaluated in detail for the LAPW basis, and the interband as 
well as the intraband contributions to the dielectric tensor are given. Results are presented for 
the metals aluminum and gold, where we crosscheck our results by sumrules. We find that the 
optical spectra can be extremely sensitive to the Brillouin zone sampling. For gold, the influence 
of relativistic effects on the dielectic function is investigated. It is shown that the scalar-relativistic 
effect is much more important than spin-orbit coupling. The interpretability of the Kohn-Sham 
eigenstates in terms of excited states is discussed. 
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I. INTRODUCTION 



Optical properties of solids are a major topic, both, in basic research as well as for indus- 
trial applications. While for the former origin and nature of different excitation processes is 
of fundamental interest, the latter can make use of them in many opto-electronic devices. 
These wide interests require experiment and theory to go hand in hand, and thus asks for 
reliable theoretical concepts. 

In solid state theory, the most successful method for the calculation of materials properties 
is density functional theory (DFT). Lattice properties like equilibrium volumes and lattice 
parameters, atomic positions, phonon frequencies and elastic constants differ from their 
experimental counterparts by a few percent only. While, however, such ground state (GS) 
properties, which are based on the calculation of the total energy, are described very reliably, 
the treatment of excited states is not rigorously justified. The main reasons are that the 
Hohenberg-Kohn^ theorem is exact only for the GS, and that the Kohn-Sham eigenstates^ 
should not be interpreted as single-electron states. Moreover, approximations have to be 
made also in the ground state calculations for describing exchange and correlation effects. 
As a consequence, the band gap problem for semiconductors is one of the most intriguing 
problems in the field. 

Nevertheless the interpretation of the KS states in terms of excited states has been suc- 
cessful for a variety of materials. Moreover, it has been claimed that the KS wavefunctions 
hardly differ from the many-body wavefunctions.—"^ It has to be pointed out that the cal- 
culation of optical properties does not go beyond the interpretation of KS eigenvalues in 
terms of the band structure. In this context, it may be hard to distinguish whether the 
shortcomings are due to the over-interpretation of GS properties or to the approximations 
used to solve the KS equations, e.g. the local density approximation (LDA). To gain more 
insight, a series of results is needed where no further approximations are introduced when 
solving the KS equations. 

To this extent we have chosen the full potential linearized augmented planewave (LAPW) 
method where no shape approximation for the potential is made. We have developed the 
formalism of treating optical properties within the random phase approximation (RPA) 
taking into account inter- as well as intraband contributions. In terms of applications we 
focus on metallic cases for the following reasons, (i) Metals don't suffer from the band gap 
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problem; (ii) metals are quite well described by the LDA, and (iii) the RPA is justified due 
to the effective screening. 

The paper is organized as follows: In Section IHl the optical response within the RPA is 
discussed, where a short comparison to time-dependent density functional theory (TDDFT) 
is made. Symmetry considerations, and the description of other optical constants obtained 
via the Kramers-Kronig relations follow hereafter. Section ITTTl is dedicated to the treatment 
of the optical response within the LAPW method. First, the basic definitions of the needed 
quantities, are given, then the momentum matrix elements are derived within the LAPW 
basis. The next section contains results obtained for the elemental metals aluminum and 
gold. Finally, conclusions are drawn on our results for the two metallic cases. The extension 
of the method with localized basis functions (local orbitals) is described in the Appendix 
together with some mathematical relations which are useful in the derivation of the formulas. 



II. THEORETICAL BACKGROUND 



A. Optical response 

The optical properties of solids are given by the response of the electron system to a 
time-dependent electromagnetic perturbation caused by the incoming light. As such, the 
calculation of these properties is reduced to the calculation of a response function that is the 
complex dielectric tensor or equivalently the polarizability. An exact expression for it is of 
course not known and one has to resort to the usual techniques of many-body perturbation 
theory to derive approximations. The first of these is the RPA. 

A particular case of the general expression for the dielectric function in the RPA is the 
well known Lindhard formula,-^ 

e(q, = 1 + ^ E ^"^^'-^-^^"^"^^'-^ , (2.1) 

where 

^(q) = (2-2) 
|q| 

is the Coulomb interaction, fie the unit cell volume, /o the Fermi distribution function, and 
£k the single particle energy. The factor 2 comes from the summation over the spin. This 
corresponds to a free electron gas in the Hartree approximation, namely, the electron-electron 



interaction is reduced to the interaction of each electron with a homogeneous self-consistent 
field. Many equivalent ways of deriving this expression have been given in the literature.i^ 
According to Hedin^ the dielectric function is defined as 

e(r, t; r', t') = 6{r - r')6{t - t') - J P{r, t; r", t')v{r" - r')dr" (2.3) 

where v is the bare Coulomb interaction and P is the polarization propagator. The simplest 
approximation for P is the RPA that corresponds to the Hartree approximation for the one 
particle Green function G^{r, t; r', t'). In this approximation, the polarization propagator P° 
is given by 

P°(r, t; r', t') = -ihG\v, t; r', t')G\r\ t'; r, t) (2.4) 

For the case of a solid with translational symmetry we can Fourier transform the equation 
for the dielectric function Eq. ()2.3|) to obtain its expression as a tensor in the reciprocal 
lattice vectors: 

ecG' (q, uj) = 5g,g' - t^(q + G)P^_g' (Q, ^) , (2-5) 

where Pq (^i{c[, uj) is the Fourier transform in the lattice of the bare polarizability defined in 
Eq. (1231). It is given by 

^^.C'(q,-) = 7^ E ^^'"'^^^ : K',Jk,q)]* M:^;(k,q) , (2.6) 

n ,n,k 1 ' -1 1 

with the matrix elements M defined as 

M^^^{k,q) = (n',k|e-'(^+^)-^|n,k + q) . (2.7) 

The expression obtained so far is the dielectric tensor that connects the total electrostatic 
potential in the solid, V, with the potential produced only by the external sources, V'^^^, 
through the expression 

VS"*(q,cu) = EeG,G'(q,^)^G'(q,^) . (2.8) 
G' 

If we denote by e the inverse of the dielectric tensor, we can invert the relation above to 
obtain an expression for the total potential 

Vg{ci,co) = J2^G,G'{ci,cu)V^^\q,u) . (2.9) 

G' 



4 



When considering perturbations produced by light, only long wavelength components are 
present in the external perturbation potential. Therefore, it is enough to assume that only 
the G = component is different form zero. In this case the response is 

VG(q,^) = eG,o(q,^)^o'"*(q,^) ■ (2-10) 

This expression shows that even in the case of an external potential with long wavelength 
variations in space, the response of the solid will have shorter wavelength components that 
are commonly called local field effects. The macroscopic dielectric constant eM, related to the 
measured optical properties in solids is the ratio between the average of the total potential 
in one unit cell, i.e. Vo(q, uj), and the external field. In this way we identify the macroscopic 
dielectric constant 

6M(q,^) = _ / , (2.11) 

by the (0, 0) element of the inverse dielectric tensor. 

The expression for the macroscopic dielectric constant given in Eq. 1)2.1111 is rather costly 
to evaluate. It involves the evaluation of the dielectric tensor with components (G, G') up 
to a certain cutoff and the subsequent inversion to obtain the (0, 0) component of the inverse 
tensor. It is a common simplification to neglect the local field effects and replace the (0, 0) 
component of the inverse by the inverse of the (0, 0) component to obtain 

e^^''\q,u) = eo,o(q,a;) = 1 - t;(q)Po°o(q, ^) • (2.12) 

We are going to neglect local field effects for the rest of the paper in order to give the 
zero order theory evaluated with the Kohn-Sham eigenvalues and eigenvectors from LAPW, 
one of the most precise electronic structure methods available. Further improvements of 
the theory can be introduced by including these local field effects or a better many body 
description of the electron system. Our results will be the starting point where no basis set 
effect can be blamed for inaccuracies. Early attempts to include local field corrections by 
Stephen Adlei^, and independently but in the same spirit by Nathan Wiser'' found that, in 
a first approximation, these effects can be described by the Lorenz-Lorentz formula with a 
renormalized polarizability. The effect turned out to be quite small for most Fermi surfaces 
studied. More recently, it has been found that local field effects are important for high 
energy excitations where localized states are involved.— 
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The Lindhard expression for the dielectric constant of Eq. ()2.ip can be obtained from the 
expression of the macroscopic dielectric constant without local field effects (Eq. ()2.12|) ) by 
assuming that the electron system is composed by only one band and the matrix elements 
are equal to one. This is the case for a simple parabolic band of free electrons. 

Since we are interested in the optical properties of solids, and the wavevector of light q 
is much smaller than any typical wavevector of electrons in the system, we need to evaluate 
the polarizability P entering in Eq. ()2.12|) in the limit of q ^ 0. The matrix elements 
involved in the expression for P given in Eq. ()2.6|) and defined in Eq. ()2.7|) can be evaluated 
for small q by perturbation theory. We will show the details in Appendix IVI B| where the 
result is given in Eq. (|6.17|) . This expression shows that the limit is different for intraband 
matrix elements, i.e., those with n' = n, and for interband matrix elements. According 
to this observation it is convenient to split the sum over n' and n into those terms with 
n' = n, corresponding to electronic transitions in the same band, and those with n' ^ n 
corresponding to interband transitions. In this manner, we can write 

eif ^ (q ^ 0, cj) = 1 + e^^'^*''^^ (q ^ 0, cj) + e^'°*'='^^ (q ^ 0, cj) , (2.13) 

where the intraband part of the dielectric constant is given by 

e{intra}(q^O^^)^_li^^![^ y /o(gn,k+q) - /o(^n,k) |MO„(k,q)|' , (2.14) 

and the interband part by 

ei-->(q-0,c.) = -lim^-i^ Y: - /"(^"■'-) |M„V(k,q)r . (2.15) 

While taking the limit of q — we can use the expression of the matrix element evaluated 
in Eq. ()6.17|) and the expansion of the band energies obtained in Eq. ()6.15|) to obtain for 
the intraband part 

n.k ^ ^ "^n.k ^ ^ 

where the derivative of the Fermi function with respect to the energy should be considered as 
a restriction to sum only over the states at the Fermi level. The momentum matrix element 
Pn,i,k is defined in Eq. ()6.13|) and evaluated for the LAPW basis set in Section fill Bl The 
interband expression is 

ei-->(q-0,.) = -lim^ VV (2.17) 



where, for a given k, c runs over the empty states and v over the occupied states. 

From these expressions, we see that the hmit is not well defined and the result depends 
on the direction of the vector q even when the limit is taken. This defines the dielectric 
constant for q ^ as a three dimensional tensor, the dielectric tensor given by 



It is clear that in the derivation of the RPA formulas presented here the unperturbed 
electronic states are described by the bare one-particle propagator G"^. In the Kohn-Sham 
approximation the many-body problem was reduced to the solution of N one-body problems. 
These single-particle problems are not representing bare particles but an effective decoupling 
through the introduction of independent quasiparticles. If the decoupling is done with the 
Hartree-Fock approximation, the scheme we just presented is conserving in the Kadanoff- 
Baynt^ sense. Here, we are going to use the RPA expression for the dielectric constant 
evaluated with the KS orbitals. This is equivalent to assume that the frequency depen- 
dent exchange and correlation kernel of TDDFT is zero, which represents a non conserving 
approximation to the problem. Nevertheless, in order to gain physical insight into the ap- 
plicability and limitations of this procedure, we will use the RPA formula utilizing the KS 
orbitals as our approximation to the dielectric constant. This procedure is also called the 
independent particle approximation because it neglects the electron-hole interaction during 
the absorption process. This should be a good approximation for metals due to the more 
effective screening of the Coulomb interaction compared to semiconductors. Comparison 
between our results and experimental data is an important step in the validation of this 
approach. It also allows to compare, when possible, with the results of more elaborated 
schemes to treat the many body effects as GW, the soultion of the Bethe-Salpeter-Equation 
(BSE), or different frequency dependent kernels in TDDFT. 

At the end of this Section, we give the expressions which will be finally computed. Most 
important is the imaginary part of the interband contribution to the dielectric tensor com- 



(2.18) 
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ponents: 



e-"*''''^(cj) = ^^^^ y Pi;n',n,kPj;n',n,k (/(^^n.k) - /(^^n'.k)) 5 (^n'.k " £n,k " t^) (2.19) 



The corresponding real parts are obtained by Kramers- Kronig transformation. 

As it is evident from Eq. ()2.18p . the intraband part is singular at u; = 0. At this point 
the plasma frequency ujpi-ij is defined by: 



2 



^ j Pi;n,nMPj;n,nM ^ (^n,k " ^f) (2.20) 



In practical calculations, a lifetime broadening F is introduced adopting a Drude-like shape 
for the intraband contribution, ej^^^^'^^uj) then reads: 

J {intra}, X pi^j ,2 21) 

Reef->(c.) = (2.22) 

Note that the imaginary part is still singular; therefore one usually works with the optical 
conductivity. The BZ integrations in Eqs. ()2.19|) and ()2.20|) are carried out by the linear 
tetrahedron method. 



B. Symmetry 



The dielectric tensor is symmetric with up to six independent components according to 
the symmetry of the crystal. Therefore the general expression for the imaginary part of e is: 



^ Im ercx Iin e^y Im e^z ^ 



Ime^y Iraeyy Im e^^ 
Ime^;^ Imey^ Im 



(2.23) 



For orthorhombic or higher symmetry only diagonal components exist. In case of cubic 
symmetry the optical properties are isotropic, i.e., there is only one independent component 
Im e^x, while for uniaxial symmetry (tetragonal, hexagonal) and orthorhombic symmetry 
there are two and three independent components, respectively. In the monoclinic case non- 
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diagonal elements occur, i.e for 7 ^ 90° the tensor takes the form: 



Im err Im e 



xy 



v 













Im 



/ 



(2.24) 



These materials are optically active or birefringent. Only in the triclinic crystal all six 
components are different. 

These symmetry considerations only hold for non-relativistic and scalar-relativistic cal- 
culations. In these cases the spin-up and spin-down elements can be evaluated separately, 
where the corresponding real parts are obtained by the Kramers-Kronig transformation (see 
below). In case of magento-optics the symmetry is reduced by the presence of a magnetic 
field as well as by the loss of time-reversal symmetry due to spin-orbit coupling. The latter 
gives rise to antisymmetric non-diagonal components, e.g. for the magnetic field parallel to 
z: 



( 



\ 





—Re e 




Ree^j^ 



xy 








(2.25) 



The corresponding imaginary parts are again obtained by Kramers-Kronig analysis. In 
case of monoclinic or lower symmetry the off-diagonal elements can therefore have two 
independent contributions, one due to a non-orthogonal crystal axis, and the other one due 
to the magneto-optical effect. 



C. Kramers-Kronig relations and optical constants 

From the imaginary part of the dielectric tensor component Im e^j the corresponding real 
part is obtained by 

, 2 /"^ cu'Imei,(a;') , , 
Re ^Sij + -p ^ ' dw'. 2.26 

Given the real part, the inverse transformation has to be used. 

With the knowledge of the complex dielectric tensor components all other frequency- 
dependent optical "constants" can be obtained. The most often used ones are the real part 
of the optical conductivity 

Re aij{u) = ^Im ey(u;), (2.27) 
9 



the loss function 

Lijicu) = -Im (-^ ] , (2.28) 
and the reflectivity at normal incidence 



, , (n-l)2 + k2 , , 

with n and k being the real in imaginary part of the complex refractive index (refractive 
index and extinction coefficient): 



, \eii{uj)\+Reeii{uj) 
nii(w) = \ (2.30) 



k,(a;) = wM^^I^5^M^ (2.31) 



The absorption coefficient is given by 

A,(u,) = MM. (2.32) 

c 

Note that the latter formulas only hold for the diagonal form of the tensor. 



D. Sumrules 

There are three sumrules which obtain information about the absorption process: 

J a {u) uj duj ^ Neff (a;') (2.33) 



J Im (^-^^ udw^ Neff {u') (2.34) 





oo 







Im ( —^1 -dLJ=- (2.35) 



The first two give an effective number of electrons contributing to the absorption process as 
a function of energy. Typically, in the low energy region the contribution to the intraband 
spectrum should sum up to the number of the outermost valence electrons. 
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III. OPTICAL RESPONSE WITHIN THE LAPW METHOD 



The LAPW basis set 



In bandstructure calculations based on density-functional theoryi the single-particle elec- 
tronic states ^nk(r) and energies e„k are described by the solutions of the Kohn-Sham (KS) 



equation* 



2m 



*nk(r) = £nk^nk(r) 



(3.1) 



with the effective potential Veff{r) being the sum of the bare Coulomb potential of the atomic 
nuclei Viatf(r), the Hartree potential Vnir) and the exchange correlation potential V^c(r)- 
In practical calculations, Eq. (jH.lll is solved via the Rayleigh-Ritz variational principle. In 
this procedure the electronic states \l'.„k(r) (KS-orbitals) are first expanded in terms of a 
physically appropriate finite set of basis functions {^k+c}. 



^nk(r) = 5^C„k(G)0k+G(r) 



(3.2) 



with G and C„k(G) denoting a reciprocal lattice vector, and the corresponding variational 
coefficient, respectively. To determine C„k(G) the ansatz ()3.2|) is inserted into Eq. ()3.ip 
followed by the minimization of the total crystal energy with respect to the variational 
coefficients. The eigenvectors and eigenvalues of the resulting matrix equation. 



y^(-^k+G,k+G' - ^nk5'k+G,k+G')C'nk(G') =0 k G BZ 



with 



k+G,k+G' = \<Pk+G 



2m 



(3.3) 



(3.4) 



-^k+G.k+G' = (0k+G |0k+G')n, (3-5) 

being Hamilton and overlap matrix, respectively, finally provide the numerical values for 
Cnk(G) and e„k- is the volume of the unit cell which in the LAPW metho d^^i^^ is 
partitioned into an interstitial region (I) and non-overlapping muffin-tin spheres (MT^) 
centered on the atomic nuclei. The corresponding basis functions are defined as, 

1 



5k+Girj 



/a 



(3.6) 
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and 

0k+G(S„ + r) = J2 K„.(k + GKir, E{) + ^^(k + G)«r(r, E{)] Y^r) (3.7) 

; ' ^ ' 

im 

= E i^O>^^.rn(f) |r| < (3.8) 

Im 

where Sq, denotes the position vector of the atomic nucleus a. The radius of the corre- 
sponding muffin-tin sphere is Ra- The product of the spherical harmonic Yi^mi^) and the 
radial function u"(r) is the solution of the Schrodinger equation for a spherical symmetric 
potential where the eigenvalue has been replaced by an appropriate energy parameter. The 
second radial function uf{r) is the derivative of the first one with respect to the energy. 
The coefficients A^(k+ G) and Bj^^(k+ G) are determined for each atom by imposing 
the requirements that the basis function has to be continuous in value and slope on the 
MT-surfaces: 



AL(k + G) = (kTG) Rl e^('^+«)«'^a,(k + G) (3.9) 

SL(k + G) = (kTG) Rl e^C^+^^^^^fo^k + G) (3.10) 

with 

ai = ][ (|k + G|i?„) ui{R^) - 31 (|k + G|i?„) u\{R^) (3.11) 

h = ji (|k + G\R^) uiiRa) - fi (|k + G\Ra) u'liRa) (3.12) 

For the explicit evaluation of the Hamilton matrix elements ()3.4j) an appropriate dual 
representation of the effective potential is needed. Within the atomic spheres the potential 
is expanded into spherical harmonics and in the interstitial it is represented by a Fourier 
series. 



B. The momentum matrix elements 

Due to the dual representation of the LAPW basis functions, the momentum matrix 
element is a sum of contributions from the atomic spheres as well as from the interstitial 
region: 

(n'k IpI nk) = ^^(n'k |p| nk)MT„ + (^'k |p| nk)i (3.13) 
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1. Contributions from the atomic spheres 



The usage of spherical harmonics in the LAPW basis set suggest to calculate the expres- 
sions (n'k \dx + idy\ nk) and (ra'k \dx — idy) \ nk) and to derive the x— and y— component 
as their linear combinations. Taking into account the ansatz (|3.2j) for the wavefunctions the 
following expressions have to be evaluated: 



k+G',k+G 



k+G',k+G 



k+G',k+G 



= (0k+G'(Sa + r) \dx + idy\ 0k+G(Sa + r)) 
= (0k+G'(Sa + r) \dx - idy\ 0k+G(Sa + r)) 
= (0k+G'(S„ + r) 0k+G(Sa + r)) 



(3.14) 
(3.15) 
(3.16) 



Expressing dx + idy, dx — idy, and dz in spherical coordinates 

dx±idy = sin^e^^<^— + -e±*^ (cos^— — 

or r \ oO sm 0(p 

oz = cosO— smt^— 

or r 09 



(3.17) 
(3.18) 



we can elaborate the matrix elements ()3.14HHTT^ . We explicitly demonstrate the procedure 
for the first component: 



a^x+iy 



k+G',k+G - \</^k+G' 



;s„ + r) 



or r \ ou sm u Oa 



0k+G(Sa + r)) 

(3.19) 



Here the first term of the operator acts on the radial coordinate only, whereas the second 
term acts on the angular coordinates only. Applying this operator to the LAPW basis 
functions ()3.8|) one obtains: 

(dx + idy) 0k+G(Sa + r) = §^^tm^''(^) sin^e^<^l^,^(f) + 



Im 



I 5^ (3.20, 

Exploiting the relations between the spherical harmonics as summarized in the Appendix, 
and omitting the arguments of the spherical harmonics for simplicity, Eq. ()3.20p becomes 



{dx + idy) 0k+G(Sa + r) = ^ 



Im 



Iw-^^^ir) + 



r 

l + l 



^Lm^l+l,rn+l + 



dr 



Wimir 



^tlYl-l^m+l ■ (3.21) 
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The X + iy component will then be 



Ra 



k+G'k+G- / 4>]i+G'{Sa + r)[dx + idy](f)i,+G{Sa + r)dr 



r 







I'm' Im 



X 



Im ' 



r) - -Vrzm(r) 
r 



it, (r) 



^l' ,l + l^m' ,m + l 



Ra 



I'm' Im ~ 



Fi^ j dnYl^,Yi,r,^^r . (3.22) 



Z.m ^ ' 



Writing explicitly the sums over /, m, /', and m', taking into account an upper limit Imax for 
/ we obtain: 

Imax I' Imax I 

ajf,x+iy 

Mr 



k+G',k+G = 5Z 5Z 5Z X] {r)Sl',l+l6m',m+l 
l'=0 m'=-l' 1=0 m=-l 
Imax I' Imax I 



J2 J2 J2 Yl Tl'm'ir)Sl>j-i6m',m+l = 
l'=0 m'=-l' 1=0 m=-l 
Imax — 1 I 

E E K^'^^'w+^S.-^-iW) (3-23) 

1=0 m=-l 

Evaluating R\^'"^~^^{r) and T/^^_-^(r) using Eqs. (j3.22|) and (j3.8p we finally obtain: 
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lmax—1 I , , , 



Z=0 m=-l 



V / 

A*^i,.n+i(k + G')S-^(k + G) u,^^{r)v!lTYdr - I f ui^,{r)ui{r)rdr) + 

R R 

G')A^^(k+ G) ( Jui+,{r)u[{rydr - I f ui+,{r)ui{r)rdr) + 

R R 

i?z+i,™+i(k + G')i^,%(k + G) f /\;+i(r)«Kr)r2rfr-/ / «;+i(r)«Kr)rdr 

Vo / 

+ G')Af+i,^_i(k + G) f / + + 2) f ui{r)ui+,{r)rdr) + 

Vo / 

A:^(k + G')Sr^i,^_i(k + G) rjuiir)u[^,irydr + (1 + 2) f uiir)ui^,ir)rdr) + 

\ R 

G')Ar+i,m-i(k + G) f/ K,(r)iiJ+i(r)r2cir + + 2) f ui{r)ui+i{r)rdr) + 

V / 

Bj-Jk + G')i?r+i,^_i(k + G) (fuiir)u[^,irydr + (/ + 2) / «;(r)«;+i(r)rrfr 

V 



(3.24) 
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with u'i{r) denoting dui{r)/dr. Analogously we derive the other two components: 

lax—l I ( 

E E 

;=o m=-l 



lmax—1 I I , , 
^k+G',k+G — )^lm 



(Ra Ra \ 

/ ui+,{r)u[{r)rHr - I J ui+,{r)ui{r)rdr + 
/ 

A*-i,^+i(k + G')S-^(k + G) (J ui^,{r)u[{r)r^dr - I J ui^,{r)ui{r)rdr\ + 

V R 

G')A^^(k+ G) ( [ ui+^{r)u[{rydr - I f ui+,{r)ui{r)rdr) + 
\ / 

i?;_"i,™+i(k + G')B-Jk + G) (fui+,ir)u[irydr - I f ui+,ir)uiir)rdr 

V / 

(Ra Ra \ 

/ ui{rX^,{rydr + (/ + 2) / ui{r)ui^,{r)rdr + 
/ 

A:5^(k + G')Sr^i,^+i(k + G) rjui{r)u[^,irydr + {l + 2) f ui{r)uUr)rdr) + 

\ R 

<J^.(k + G')^r+i,m+i(k + G) (/ K,(r)K;+i(r)r2cir + (Z + 2) f ui{r)ui+i{r)rdr) + 

V / 

^i(k + G')i?r+i,^+i(k + G) (fuiir)u[^,irydr{l + 2) / «;(r)«;+i(r)rrfr 

\ 



(3.25) 



Imax—l I J 

"^k+G',k+G — E E ] -^Zm 

1=0 m=-l ^ 



A+i,m+i(k+ G')A^„(k+ G) M,+i(r)KKr)r2cir - / f Ul+^{r)ul{r)rdr) + 

\o / 

A*Ti,m+i(k + G')i?^^(k + G) (f ui+^ir)u\irydr - I f ui+^{r)uiir)rdr) + 

\0 / 

^;^i,n.+i(k+ GOA^^(k + G) rjui+^{r)u[irydr - I f ui+,{r)uiir)rdr) + 

\ / 

^m,m+i(k + G')S^^^(k + G) rfui+i{r)u[{rydr - I f ui+^{r)ui{r)rdr 

V / 

, ^(6) 

A*,^(k+ G')^r+i,m+i(k+ G) (fui{rX^,{rydr +{l + 2) f ui{r)ui+i{r)rdr) + 

A;^(k + G')5r+wi(k+G) (fuiir)u\^,irydr + {l + 2) J uiir)ui+^ir)rdr) + 

\ / 

B:^ik + G')Ar+i,.n+i(k + G) (fuiirX^.irydr +(1 + 2) / ^KO^w(r)rdr) + 

R i? 

K^i^ + G')^r+i,m+i(k + G)( [ ui{r)u\^,{rydr + + 2) [ ui{r)ui+,{r)rdr 

\ 

(3.26) 
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To obtain the atomic-sphere contributions to the momentum matrix elements, Eqs. ()3.14t - 
I3.16|) have to be multiphed with the variational coefficients and summed up over all basis 
functions: 



x+iy j_a ^x-iy \^ 

k+G',k+G ' ^k+G',k+G;' '^nky^) 



G'G 



in 



'k \Vy\ nk)MT. = ^ 5^ C:,^{G') {-KV&^^^a *k+^',k+G) C^nk(G) 



G'G 



(n'k |V,| nk)MT. = Yl C:.k(G')"$^+G',k+GCnk(G) 



(3.27) 



G'G 



2. Contributions from the interstitial region 



With the planewave basis ()3.6|) the interstitial contribution to the matrix elements ()3.13p 
can be easily worked out: 



{n' 



'k|Vhk)i = ^J2 Q'k(G')C„u(G) [ e^(«'-«)''cir 



(3.28) 



The integration in Eq. ()3.28|) over the interstitial region is carried out by integrating over 
the whole unit cell and subtracting the integral over the atomic spheres (see Appendix). 
This procedure finally leads to 



(n'k I V| nk)i = ^ J] (k + G) C„k(G) 



c:,k(G) fie-E^- -E^^" 



E ^"'k(G 
G'^G 



,^ Sin (|G' - G| R^) - (|G' - G| R^) cos (|G' - G| R^) ,(g-g')s„ 



(|G'-G|i?„) 
where Va is the volume of the atomic sphere a. 



(3.29) 



IV. RESULTS 
A. Aluminum 

As the first example fee aluminum has been chosen. It has been intensively studied in 
Iiterature,'^^*^^»^^*i2i21i where the most puzzling experimental problem at that time was the 
distinction between interband and intraband contributions to the absorption process. Here, 
it serves as a test case mainly in terms of convergence with respect to the most important 
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FIG. 1: Imaginary part of the frequency dependent dielectric function of Al for a series of k point 
meshes. In the inset the dependence of the plasma frequency (in eV) on the k point sampling is 
depicted. 

convergence parameters. While the results did not turn out to be sensitive to the number 
of LAPW's, there is a crucial dependence on the BZ sampling as already found by Lee and 
Chang.— Fig. ^ shows the interband contribution to the imaginary part of the dielectric 
function (Eq. I2.19|) for a series of k point meshes. The peaks do not only get sharper with 
increasing number of k points, but they also exhibit a pronounced energy shift. Only with 
the most dense mesh quantitative agreement with experimental data is achieved.— The inset 
of Fig. n exhibits the plasma frequency (Eq. I2.2()j) as a function of the number of k points 
in the irreducible part of the Brillouin zone (IBZ). Also in this case highest accuracy is 
needed to reach the converged value of Upi = 12.6 eV. This sensitivity can be understood 
by the high symmetry of the crystal structure and the simplicity of the material. Since 
the main contributions stem from certain regions of the BZ/^ a refinement of the mesh 
there can dramatically change this contribution and - due to the high weight in the BZ - 
the total spectrum. In contrast, in more complex materials the optical absorption usually 
represents a sum of many different contributions from several regions in the BZ and many 
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band combinations. Therefore a change of one or some of these contributions cannot affect 
the total spectrum that much. Fig. |21 illustrates the separation of interband and intraband 
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FIG. 2: Loss function for Al (full line) calculated with 4735 k points in the IBZ. The dashed 
(dotted) line shows the loss function of the pure intraband (interband) contribution. 

contributions to the optical spectra. For this purpose the total loss function is plotted (full 
line), but also the corresponding functions if there was interband (dotted line) or intraband 
(dashed line) contributions only. Judging from the shape of the curves all of them could be 
interpreted as free electron behavior but with very different plasma frequencies. Therefore 
the plasma frequency was for a long time thought to be 15.2 eV which is the peak position 
of the total spectrum. For a more detailed discussion of the problem see Ref. Q. Finally 
the sumrules provide another test for the quality of the calculations. To this extent, all 
three sumrules as given in Eqs. ()2.33l - 0.35|) have been tested. In Fig. El the results of the 
corresponding integrations are depicted as a function of energy, i.e. the upper integration 
limit. The theoretical results for sumrules 12.331 and 12.341 excellentlv reproduce the region 
where experimental data are available, while the third curve approaches 7r/2 within 0.1 %. 
Having a closer look to sumrule (1), the kink in the low energy range indicates the setting 
in of the interband contributions. 
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FIG. 3: Sumrule check for the optical spectra of Al by carrying out the integrations given in Eqs. 
12.331 fl). 12.341 (2). and 12.351 (3) as a function of the upper integration hmit (energy Tluj in eV). The 
calculated values are given for two different numbers of k points in the BZ, i.e. 165 (dashed lines) 
and 4735 (full lines). The corresponding experimental data (open circles) are taken from Ref. Ilfil . 
The dotted line indicates 7r/2 for comparison with the third sumrule (3). 

B. Gold 

The electronic bands have been studied very intensively in the liter at urej-^^'^SiSSiSi and 
much attention has been paid to relativistic effects on the band structure and henceforth 
the optical properties.— But the optical spectra have not been investigated in detail from 
the first-principles side. Therefore they shall be studied in more detail here. 

A lattice parameter of 7.66 a.u. has been taken which is the optimized value with respect 
to a scalar-relativistic LDA calculation. The self-consistency cycles have been performed 
with 65 k points in the IBZ. The BZ integrations for obtaining the density of states (DOS), 
joint DOS, and the optical properties used in the figures were carried out with 1240 points 
in the IBZ. (For convergence tests see below.) In order to investigate relativistic effects we 
have performed non-relativistic, scalar-relativistic, and relativistic calculations, where in the 
latter spin-orbit coupling has been treated by a second variational scheme. In all cases the 
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core states are calculated fully relativistically by solving the Dirac equation. We treat the 
semi-core states 5s, 5p and 4/ in the valence region by local orbitals. Therefore these low 
lying states cannot be expected to be very well described in the non-relativistic case. Our 
procedure, however, has the advantage that the optical properties can be handled by the 
formalism described above up to high excitation energies which is needed for the Kramers- 
Kronig analysis as well as for the test of the sumrules. At the same time the influence of 
relativistic effects can be studied in detail. 

Also in the case of gold, we have paid attention to the convergence of the optical spectra 
with respect to the k point sampling. It turned out that the density of the k point mesh is 
much less important than in aluminum. Already rather course meshes give all results with 
high accuracy. For example, the plasma frequency of the relativistic calculation changes by 
0.1 eV only, i.e. from 8.966 to 8.867, when increasing the number of k points in the IBZ 
from 165 to 11480. In fact, already with 47 k points a value of 8.866 eV is obtained. For the 
interband contributions, the spectra above 2 eV can hardly be distinguished when plotted 
for different BZ samplings. Only the first interband peak (below 2 eV) in the real part of the 
dielectric function is somewhat broader with more spectral weight at lower energies when the 
number of irreducible k points is reduced below 1240. In Fig. |3]the imaginary part of the 
dielectric function arising from interband transitions is shown for the three cases mentioned 
above. While there is hardly any difference between the scalar-relativistic calculation and 
the treatment including spin-orbit coupling the biggest change is observed when taking into 
account the scalar-relativistic effects which are due to the Darwin shift and the mass-velocity 
term. 

In order to analyze the origin of theses effects the corresponding band structures are 
shown in Fig. One striking scalar-relativistic effect is the broadening and upward shift of 
the filled 5d bands, which is much bigger than the splitting of some states due to spin-orbit 
coupling. Also the 5s states lying at -9 eV at the F point in the non-relativistic case exhibit 
an increased band width by more than a factor of two. Relativistic effects of the 6s and 5p 
states in the conduction band are highlighted in the figure by indicating the predominant 
band character along the high symmetry lines X-W and W-K. The 5/ states dominate the 
band structure above 12 eV in the non-relativistic case, and these bands are pushed up by 
roughly 1.5 eV. All these effects can be seen more clearly in the total and site-projected 
densities of states depicted in Fig. IHl 
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FIG. 4: Interband contribution to the dielectric function of gold: real (grey lines) and imaginary 
part (black lines) of Im e from non-relativistic (top), scalar-relativistic (middle) and relativistic 
calculations (bottom). 



The first half of the spectrum is dominated by transitions between d and p like states 
which generally are shifted to lower energies due to the upward shift of the d band when 
relativistic terms are included in the calculations. In particular, the non-relativistic Im e 
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FIG. 5: Relativistic effects in tlie band structure of gold: non-relativistic (rigfit), scalar-relativistic 
(middle) and relativistic calculations (left) have been performed. The main Z-like band character is 
highlighted by symbols as indicated in the legend. The thin lines without symbols represent states 
with more than two different significant contributions. 

centered around 5 eV is created by 5d to 6p transitions with the partially filled band of 
mainly p character providing the final states. This peak is shifted down by 2 eV in the 
scalar-relativistic case due to the nearly rigid upward shift of the 5d bands. The increase 
of the spectral weight can be traced back to the l/o;^ behavior of Im e. The double peak 
structure around 6-8 eV visible in the relativistic cases come from transitions from the three 
highest lying valence levels to the first two unoccupied bands. They also exhibit mainly d to 
p character, but due to some p admixture in the high lying valence states and a considerable 
amount of d character in the final states also finite p to d contributions are present. Also 
the d to f transitions around 20 eV are shifted to lower energies because the / states move 
up only half as much as the d states. Including spin orbit coupling only very minor changes 
are obtained. In particular, the d to f transitions are broadened reflecting a splitting of 
the d bands as can be seen in the band structure and the symmetry projected DOS. The 
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Energy [eV] 

FIG. 6: Total and site-dependent densities of states of Au from a non-relativistic (dark grey), a 
scalar-relativistic (light grey) and a relativistic calculation (black), spin-orbit coupling 

transitions in this case exhibit an asymmetry with respect to the Fermi leveL In all cases 
the initial states are mainly the d bands between -10 and -1 eV, while the final states range 
from the Fermi level up to very high energies. As already pointed out by Christensen24 the 
matrix element effect is relatively small. The imaginary part of the dielectric function is 
already very well estimated by the joint density of states as can be seen in Fig. Since 
the dielectric function suppresses higher energy transitions due to its l/cj^ behavior (Eq. 
I2.5|l we have divided the joint DOS by the energy squared. Only a few features there are 
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FIG. 7: Joint density of states of gold divided by energy squared. The data are shown for the 
non-relativistic case (dark gray full line), the scalar-relativistic case (light grey full lines) and the 
relativistic case (black scattered line). 

canceled by selection rules, like the sharp peak at 4.5 eV arising from the parallel bands 
of d character at -2 eV and s character at 2.5 eV as clearly seen along the line W-L in 
the BZ. In the non-relativistic case up to 13 eV, the magnitude is much too small and the 
peaks are found at too high energies as already seen in the optical spectra. Comparing 
our results to experimental data (Fig. |SJ, the overall agreement is good. There is two main 
points which have to be considered. First, the experimental spectrum is in general smaller in 
magnitude. This seems to be a common problem which can most probably be traced back to 
a diminished reflectivity at the sample surface.-^^^ Second, although all experimental features 
are reproduced by theory, the peak positions appear to be slightly too low in the latter. 
Since most of the transitions contain d states either as initial or as final state, the width of 
the d band is presumably the source of this shortcoming. The number of effective charges 
Nf,ff probed by the two sumrules (Eqs. ()2.33p and ()2.34|) behave quite differently as can be 
seen in Fig. 121 While the integration of Ime reflects the different magnitudes of theoretical 
and experimental data, the integration over the loss function is much less sensitive showing 
very good agreement between theory and experiment in the whole energy range. 
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FIG. 8: Imaginary part of the dielectric function for gold including a life-time broadening of O.leV 
obtained by a relativistic calculation compared to the corresponding experiment taken from Cooper 
et M 

Before the over-interpretation of KS states are blamed for the small discrepancies, other 
reasons should be investigated. From the theoretical side they can be found in the correlation 
effects being underestimated by the LDA as discussed above, but also in the rather poor 
basis set of linearized methods for the conduction bands. Furthermore, for the high energy 
photons the g = limit is not justified any more and thus a finite momentum transfer should 
be considered in future work. But also the experimental side gives rise to shortcomings, as 
e.g. the Kramers-Kronig analysis requires a continuation of the reflectance spectrum. We 
need to include Im e up to more than 60 eV in the KK transformation in order to obtain 
Re e and henceforth the loss function up to 30 eV reliably. 

Although the changes of the band crossing the Fermi level are hardly visible from the 
band structure, the plasma frequency changes from 9.02 to 9.51 when taking into account 
scalar-relativistic effects and becomes 8.86 when spin-orbit coupling is introduced. This 
sensitivity of this quantity is a measure for the reliability of the Fermi surface by the fuUy- 
relativistic calculation, ujpi excellently agrees with its experimental counterpart of 8.83.— 
An explanation can also be given for the fact that the plasma frequency hardly depends on 
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FIG. 9: Sumrule test for the optical spectra of gold obtained by a relativistic calculation compared 
to the corresponding experiment taken from Cooper et al.^^ The numbers in brackets refer to the 
corresponding sumrules (Eqs. (|2.33l and l2.34() . 

the k point mesh. Since the only band contributing to the effective number of free carriers 
has a linear shape around the crossing with the Fermi level the linear tetrahedron method 
used in our work turns out to be a very effective tool for the Fermi-surface integration. 



V. CONCLUSIONS 

We have provided a scheme for the calculation of linear optical properties of solids in the 
random-phase approximation. It is worked out within the LAPW method which does not 
require approximations in the solution of the KS equation and thus is a good starting point 
for the evaluation of the KS states in terms of band structure and excited states. 

Aluminum, which is a prototype for a free-electron-like metal, exhibits excellent agree- 
ment with experiments in terms of plasma frequency, peak positions of the interband tran- 
sitions, as well as sumrules. This agreement can, however, only be achieved when the BZ 
integration is carried out with extreme accuracy. Thus highest precision is in general needed 
to draw conclusions about the quality of the KS states. 
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For gold, the optical properties are less sensitive to the BZ sampling. The plasma fre- 
quency is in excellent agreement with experiment when relativistic effects are taken into 
account including spin-orbit coupling. The latter is, however, less important for the in- 
terband transitions than the scalar-relativistic contributions, i.e. the Darwin shift and the 
mass-velocity term. Shortcomings can at least be partly traced back to the approximations 
used for the ground state calculation, i.e. the LDA which does not properly account for 
correlation effects and thus cannot bring out the d band width and location well enough. 
Other possible weak points are the insufficient basis set for the unoccupied states, the g = 
limit for high energy transfer, as well as the missing local-field effects. The latter, however, 
turned out to be of minor importance for the calculation of hot-electron lifetimes.— All 
these facts have to be accounted for before a final statement about the interpretability of 
KS states for a certain material can be made. Nevertheless, in both cases studied here, they 
provide at least a very good first approximation to single state and excited state energies, 
respectively. 
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VI. APPENDIX 

A. k.p Perturbation Theory 

In the position representation, the one-electron wavefunction (r | n, k) = "^n.kli") is the 
solution of the Schrodinger equation 



Since the crystal potential V^(r) has the periodicity of the lattice, using Bloch's theorem we 
can write 




+ V{y) V^„,k(r) = en,kV^n,k 



(6.1) 



(6.2) 
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The periodic part of the wavefunction u„,k(r) obeys 



where 



(P + 
2m 



-rV{v) 



(6.3) 



(6.4) 



We want to obtain an expression for Un,\L+q{v) and £n,k+q vahd for small q in terms of the 
values for q = by perturbation theory. The equation for w„^k+q(r) is 



-^^k+q 1in,k+q(r) — £n,k+q 1in,k+q(r) , 



with 



2^2 



k+q 



/ip ■ q /i k ■ q h q 



+ 



+ 



m m 2m ' 

where the last three terms can be treated as a perturbation. 
At this point it will be useful to introduce the notation 



(6.5) 



(6.6) 



and 



n', k' 



O 



(r I n, k) = Un,k{r) 



(6.7) 



(6.8) 



for any operator O with the integral running over the unit cell of volume fl. 

According to this notation, the wavefunction to first order in q for a non-degenerate state 

is 

I 1 , X I 1 N , I / 1 X (^'>k|i^k+q-^k|n,k) 
In, k + q) = In, k) + > In ,k) ^ (6.9) 



^n.k ^n'.k 



Only one term from the perturbation Hamiltonian gives a contribution different from zero 
to the first order correction of the wavefunction 



,k + q) = |n,k) + J]|n',k).^ ' ''"^ ^1 ^ 



^n.k £-n',k 



(6.10) 



To linear order in q the expression for the energy is 



hp ■ q k • q 



n,k . 



^n,k+q — ^n,k + I n, k 

* m m 

Since the states are normalized this expression can be written as 

£n,k+q = £n,k + — [(^i, k |p| n, k) + k] • q 



(6.11) 



(6.12) 
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The momentum matrix element evaluated in Section IIIIBI can be expressed as 

P/,n,k = {I, k IpI n, k) = 6i^n hk+ {l,k |p| n, k) . (6.13) 
In terms of this definition, the wavefunctions and energies to first order in q are given by 

\n, k + q) = |n, k) + - V \n', k) P"''"'*^ ■ q (6.14) 

and 

h 

£n,k+Q = £n,k -\ Pn.,n.k (6.15) 

m 

B. Matrix elements for small q 

In this appendix we use the k-p expressions developed in Appendix IVI Al to evaluate the 
matrix elements for q ^ 0. We first write the matrix elements in terms of the periodic part 
of the wavefunction 

M°„(k,q) = (/,k|e-*'i"|n,k + q) = (/,k| n,k + q) . (6.16) 

Using the expression of the wavefunction for small q of Eg. ()6.14j) . the matrix element at 
lowest order is 

M°„(k, q ^ 0) = + (1 - 5,,„)- -P^^ ■ q (6.17) 

C. Relations between spherical harmonics 

Below, several well known relations between the spherical harmonics are given, which are 
needed to evaluate the momentum matrix elements: 

e+^^sin^y,,^ = + i'S^^-i.^+i (6-18) 

e-'^ sin OYi,^ = F^^Yi+i^m-i + F.^^l^-i,™-! (6.19) 
cos^r^,„ = + i^S>^-i,m (6.20) 
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with 



(1) _ _ / (/ + m + l)(/ + m + 2) 
(2/ + l)(2/ + 3) 



.(2) _ / (/-m)(/-m-l) 
(2/ - 1)(2Z + 1) 



Lm 



p(6) 



'(/ - m + -m + 2) 
(2/ + l)(2/ + 3) 

'(/ + m){l + m - 1) 
(2/ - 1)(2/ + 1) 

'(/-m + l)(/ + m + l) 
(2/ + l)(2/ + 3) 

(/ — m)(/ + m) 
(2/-l)(2/ + l) 



6.21) 
6.22) 
6.23) 
6.24) 
6.25) 
6.26) 

6.27) 
6.28) 
6.29) 



D. The step function 

The analytic evaluation of a planewave integral over the interstitial region as needed 
in Eq. ()3.28|) is carried out by integrating over the whole unit cell and subtracting the 
contributions of the atomic spheres, which are 

[ gi(G'-G)r^j. _ J (Q 30) 

/ I OT/ sin(fl„|G'-G|)-cos(ii^|G'-G|) j(G'-G)S^ P' ^ P V • y 

MT. I '^^^ {Rc.\G'-G\f ^ ^ 

It utilizes the Rayleigh-expansion of a planewave in terms of spherical harmonics: 

e^Gr ^4^g^GS.^^«^.^^|r_s,|G)F4(G)Fz™(r^„) . (6.31) 

lm 
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E. Local orbitals 



The extension of the LAPW basis set by locahzed orbitals^^ was introduced in order to 
describe semi-core states, i.e. those low-lying states which reach out of the atomic sphere, 
on the same footing as valence states. The corresponding basis functions reads inside the 
atomic sphere 

0k?„.(S. + r) = 

[AL(k + GL°)<(r, El) + Br^ik + GL^)«r(r, E^) + Cr^(k + GL°)<(r, E^)] F,,^(f ) 

(6.32) 

and is zero outside. The coefficients Af^ik + Gf^), B^^ik + Gf°), and Cj'^ik + Gf^) are 
determined by choosing 0^,?™ spatial derivative to vanish at the sphere boundary, 

and by the normalization condition for 0k?m- The additional trial energy E[''-' corresponds 
to the energy of the semi-core state, and for each localized basis function one specific Gf"^ 
is chosen.— The wavefunction can then be written as 

*nk(r) = J2 ^nk(G)0k+G(r) + J2 CZmiGi^)<PtU'^) (6.33) 

G Im 

Analogously to Eqs. (|3.14flTTH|) matrix elements between LAPW's and LO's 

"^k^^k+G = (0k?r„(Sa + r) \dx + zdy\ 0k+G(S. + r)) (6.34) 
"<^krJk+G = (0kt(S. + r) \dx - idy\ 0k+G(S. + r)) (6.35) 
"'^'L^.k+G = «?™(S. + r) \dz\ 0k+G(S. + r)) (6.36) 

and between LO's and LO's have to be defined: 

(0k?m'(S. + r) \dx + idy\ 0^,1(S, + r)) (6.37) 
(0k?m'(S« + r) \dx - idy\ cl>i%{S^ + r)) (6.38) 
(C™'(S„ + r)|9z|0^,l(S„ + r)) (6.39) 



\d'm' .\dm 
\d'm' .\dm 



kZ'm' kZm 
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The atomic sphere parts of the momentum matrix elements given in Eq. ()3.27|) have to be 
supplemented by contributions from the local orbitals: 

G',G 
lm,G 

I 1 ^*LO (r^LO \ (a^x+iy .a ^x-iy \ ^LO /(~tLO\ 

"T 2 ^n'k,l'm'\^l'm' ) [ ^kl'm' ,]ilm ~^ ^kl'm' ,klm) ^nk,lm^^lm ) 

I'm' ,lm 

(6.40) 

{n'k I V,| nk)MT. = ^ C:,k(G') {''KV&'M+a ^k+cl'.k+c) C^nk(G) 

G',G 

I ^*LO (r^LO\ fa^x+iy a ^x-iy \ ^ 

+ ^ '-'n'k,/ml*^/m J V ^kim,k+G ~ ^kim.k+G^ '-^nkl'^J 

lm,G 

_l_ ^ ^*LO (CLO \ (a^x+iy a ^x-iy \ ^LO ((~^LO\ 

+ ^ 2-^ ^n']i,l'm'V^l'm' ) \ ^ki'm',kim ^kZ'm',k«my' '-^nk.Zml'^Zm J 

(6.41) 

(n'k I V,| nk)MT. = E G:,k(G') "$^+G',k+G C^nk(G) 

G',G 

+ ^ ^nn^lmi^f'm) "^k/m,k+G C'„k(G) 

lm,G 

+ CrJk!l'm'i^l'm') "^L'm'.kZm ^nkM^^lm) (6.42) 

The evaluation of Eqs. (j6. 34116. 35|1 is straightforward analogous to Eqs. (j3. 24113. 20|) . 
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